The Annals of Applied Statistics

In biomedical research, multiple types of high-dimensional data, such as genomic, transcriptomic, proteomic, and metabolomic data, are increasingly collected from the same subjects. Integrating these multiple data views can improve prediction by exploiting shared or complementary information across the views. Cooperative learning provides an agreement-based framework for multiview supervised learning by encouraging predictions obtained from individual views to be similar. However, the original framework assumes independent observations and therefore does not account for clustered structures, such as repeated measurements obtained from the same subject. To address this limitation, we propose Cooperative Learning with a penalized Linear Mixed Model (CL-pLMM) for high-dimensional multiview data with a clustered structure. CL-pLMM replaces the ordinary prediction loss in cooperative learning with a covariance-weighted loss that accounts for within-cluster dependence, while retaining the agreement penalty between views and a Lasso penalty for variable selection. We further show that its objective function can be represented as a penalized linear mixed-effects model applied to augmented data, allowing existing estimation procedures to be used. The performance of CL-pLMM is evaluated through simulation studies under various signal and dependence settings and an application to longitudinal proteomic and metabolomic data for predicting the time to spontaneous labor.

Selecting stochastic gene-expression models from single-cell counts requires balancing goodness of fit against unnecessary mechanistic complexity. The probability-generating-function-based Bayesian information criterion (PGF-BIC) combines covariance-weighted fitting in generating-function space with a complexity penalty, allowing candidate models to be compared without reconstructing their full count distributions. However, its conventional zero-threshold rule does not account for sampling uncertainty in the fitted score difference and may therefore favor overly complex models in finite samples. To address this limitation, we develop an uncertainty-aware PGF-BIC rule that selects the more complex model only when its score advantage exceeds a data-driven threshold. We use Cantelli's one-sided inequality to motivate a selection margin expressed in terms of a standard deviation. To determine this scale, we use influence functions to quantify sensitivity to small perturbations in the data distribution and obtain a first-order description of sampling fluctuations. The resulting variance estimate accounts for variability in both the empirical probability generating function and the estimated covariance weights, yielding a data-driven threshold for assessing the complex model's score advantage. A Poisson versus Bursty benchmark shows that the calibrated rule reduces incorrect selection of the more complex model. The calibration requires neither resampling nor additional optimization, incorporating sampling uncertainty into model selection while retaining the computational efficiency of PGF-BIC.

Data augmentation is now a standard device across capture--recapture and occupancy analysis: adding a fixed number M of all-zero encounter histories replaces a model of unknown dimension with one of fixed dimension. Although M is often treated as a computational tuning choice, it also specifies a finite superpopulation and hence a binomial support constraint on the number of undetected individuals. In a Bayesian implementation that constraint appears as an induced prior; in a likelihood implementation it is the same finite-support assumption reached by another route. That the Bernoulli specification for the inclusion indicators induces a binomial prior on abundance is established (Schofield & Barker 2014); our concern is what that choice costs in estimated uncertainty. We develop the argument using a simple closed-population abundance estimation problem. We show that augmented occupancy and Huggins conditional-likelihood analyses give numerically identical point estimates of N once M is sufficiently large. Their uncertainty estimates, however, need not agree. We distinguish two sources of discrepancy. First, when M is small relative to the number of undetected individuals, the finite binomial ceiling truncates the likelihood or posterior and suppresses uncertainty. Second, once that ceiling no longer binds, Taylor-series (Delta-method) approximations still understate variance, because the quantity of interest is a strongly non-linear function of the estimated parameters and local linearization does not reproduce its curvature. Gauss-Hermite quadrature on the unconstrained logit scale recovers much of the shortfall and approaches the MCMC posterior benchmark, though a small residual remains that does not close as M grows, reflecting the distinction between asymptotic likelihood theory and finite-sample Bayesian inference. Neither mechanism is peculiar to abundance estimation: the first follows from the augmented representation itself, the second from any derived quantity that is a non-linear function of estimated parameters. We close with framework-specific guidance for choosing M.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Estimating the correlation between subject-specific rates of change is a common objective in longitudinal studies. A natural approach is to fit ordinary least squares regressions per subject and correlate the resulting slope estimates directly. Under complete data, this slopes-first estimator is algebraically equivalent to a contrast-based plug-in estimator derived from the joint variance-covariance matrix of the two signals across time. Under independent dropout, however, the slopes-first approach implicitly applies heterogeneous weights across missing-data patterns, introducing bias whenever the cross-covariance between the two signals is not constant over time. We formalise this result and propose an unbiased alternative that estimates the joint covariance matrix before applying the contrast. We further extend our results to include covariates as well. Simulations confirm that the slopes-first estimator is biased under non-constant cross-covariance while the proposed alternative remains unbiased. Further, we show in a real-world example that the two estimators can lead to different scientific conclusions and recommend that researchers working with longitudinal data inspect the missing data pattern to make an informed decision about the estimator of choice.

Spatial expansion is a common feature of many ecological processes, including biological invasions, range shifts, population recovery, and the spread of pathogens and human activities. Although such expansion often arises from diffusive processes, there remain few general statistical methods to quantify these dynamics. Here, we present a diffusion-enhanced spatio-temporal model (DESTM) that approximates diffusive process using Gaussian Markov random fields (GMRF). The model works by incorporating the statistical interaction between spatial diffusion and temporal lags into a joint precision structure of the latent field. To demonstrate, we simulate data from a diffusion movement process using a continuous-time Markov chain (CTMC), and then approximate this mechanistic process using the DESTM. We show that separable spatio-temporal models can lead to large bias in estimates of the local density. In contrast, the DESTM did not appear to bias inferences, and reverts to the separable model when diffusion is not apparent. Finally, we apply the model to two large-scale examples of spatial expansion: (1) the historical spread of the Japanese distant-water longline fleet across the Pacific, representing anthropogenic spatial expansion in the ocean, and (2) the recovery and geographic expansion of Bald Eagle populations across North America. Estimates from both case studies, supported by model comparisons and ecological plausibility, suggest that the DESTM provides a better description of the underlying expansion dynamics, although improvements in forecast performance were more apparent for the Japanese longline fleet. We highlight that the model provides a fast and efficient framework for modeling spatial expansion dynamics, with the potential to represent increasingly complex ecological dynamics through future extensions. We propose that future research explore this framework as a quantitative tool to explain how spatial expansion occurs.

Advances in single-photon sensitive detectors are rapidly expanding the adoption of photon-counting fluorescence microscopy. Under the Poisson photon-counting statistics, iterative Richardson-Lucy (RL) -type algorithms are statistically optimal for image deconvolution but suffers from a fundamental semi-convergent behaviour: prolonged iterations inevitably amplify noise, requiring heuristic early stopping or regularisation typically based on assumptions about object morphology. Here we present a morphology-agnostic regularisation framework for RL deconvolution that exploits independent noise realisations instead of structural object priors. Preserving the Poisson statistics of the acquisition process, we formulate the Regularized by Noise (RbN): a regularized variational framework and its corresponding iterative minimization algorithm that exploit the statistical consistency of independent noise realisations to distinguish reproducible image features from stochastic noise. The regularisation strength is selected automatically using the Poisson residual whiteness principle, resulting in a fully data-driven reconstruction without heuristic parameter tuning. Photon-timing-resolved systems naturally provide the independent noise realisations exploited by our framework, whereas computational photon splitting provides a statistically equivalent implementation for photon-counting systems. We validate the approach experimentally using photon-timing-resolved confocal microscopy and image scanning microscopy across eight morphologically distinct subcellular targets, and further demonstrate its applicability to photon-counting microscopes through computational photon splitting. Across imaging modalities, detector technologies, biological structures and signal-to-noise regimes, our method eliminates RL semi-convergence, removes sensitivity to the stopping criterion, and consistently outperforms conventional RL while preserving fine structural detail. More broadly, our results establish independent noise realisations as a general source of morphology-agnostic regularisation. Although demonstrated here for SPAD-based laser-scanning microscopy and Poisson statistics, the underlying principle could be extended to other imaging modalities and, more generally, to statistical inverse problems through appropriate noise-specific formulations.

Multimodal Gaussian process factor analysis provides a flexible framework for dimensionality reduction in temporally or spatially resolved omics data. Existing approaches, however, typically rely on pre-specified Gaussian process kernel families and do not explicitly separate each latent factor into a component capturing gradual, smooth variation and a complementary component capturing fine-scale, non-smooth variation. Here, we present MOFTy, a Bayesian multimodal factor analysis framework based on numerical information field theory (NIFTy) that replaces fixed kernel families with the flexible correlated field model in NIFTy and enables explicit additive component separation within each latent factor with quantified uncertainty. NIFTy has been successfully applied to high-resolution Bayesian imaging in astrophysics and facilitates scalable, curvature-aware variational inference for efficient posterior approximations. We validate MOFTy on simulated data; applications to published multi-omics data demonstrate that MOFTy disentangles latent spatial structures by separating smooth gradients from localized fine-scale heterogeneity in human glioblastoma and recovers cross-modal patterns in a mouse gastrulation dataset.

Estimating differential gene expression is a common task in RNA-seq. Current methods mostly rely on normalization to reduce variance and increase accuracy. These methods are widely used and provide invaluable information about transcriptomic changes between biological conditions. However, widely used normalization methods are known to distort estimates of transcript differential expression when a majority of genes are upregulated or downregulated, or when the total RNA content per cell changes. Despite the presence of global expression shifts in a number of contexts, few methods exist that can provide accurate normalization and estimate linear models under this context without spike-in controls. Here, we present \textbf{GLORB} (\textbf{G}\textbf{L}\textbf{O}bal-shift \textbf{R}obust \textbf{B}ayesian model), a method for estimating generalized linear models under global upregulation. We develop two models that are able to recover differentially expressed genes and linear model coefficients with lower distortion of results. We show that our model's method of accounting for library size variance is consistent with DESeq2's median of ratios, and edgeR's trimmed mean of M-values under conditions when a minority of genes are upregulated and outperforms them under circumstances when most genes are either increased or decreased between groups. Finally, because our method does not rely on calculating geometric means for each gene it is able to work in datasets with much higher sparsity.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

Somatic mutations in sequenced tumors are often detected with imperfect sensitivity. They can be missed because tumor heterogeneity results in some mutations being absent in a sample at hand or present at very low frequencies which, combined with technical limitations like insufficient sequencing depth, reduces their detectability. It becomes apparent with multi-regional sequencing that some mutations confined to specific areas can be overlooked in single-region analyses. Adjusting for such sensitivity becomes critical when mutation frequencies between groups of subjects are compared, such as responders and non-responders to a given therapy. Otherwise, the estimate of the mutation prevalence, or of the difference in prevalences between the groups, can be biased, leading to inflated type I error. Here, we propose a weighted approach to adjust for imperfect sensitivity of the mutation detection and compared it to the unweighted approach. We assume that a surrogate of mutation sensitivity is observed only among patients with called mutation, and we derive its unconditional mean. Such adjustment leads to less biased estimates of mutational frequencies and better control of the type 1 error.

Natural killer (NK) cells display striking variability in their cytotoxic responses, but whether this variation reflects random events, stable differences between cells, or effects of previous encounters remains unclear. We developed a Bayesian framework that uses single-cell interaction histories to disentangle these sources of variability. The framework models target encounters and killing decisions to quantify population heterogeneity, determine the sample sizes needed to distinguish competing mechanisms, and separate stable cell-to-cell differences from history-dependent behaviour. Synthetic data established when these mechanisms can be reliably distinguished in practice. We tested the framework using time-lapse imaging of NK cells exposed to rituximab or CC-96673, an antibody co-targeting CD20 and CD47. Despite similar overall killing, the framework identified distinct underlying responses: rituximab mainly increased mean killing rate, whereas CC-96673 reduced cell-to-cell variation. Event-count analysis supported continuous population heterogeneity, while ordered interaction histories revealed that both stable differences in killing propensity and previous encounters shape cytotoxic decisions. To make the framework directly usable, we introduce BARRACUDA, an open-source web platform and Python package implementing these analyses, including donor-aware extensions, for single-cell cytotoxicity datasets.

Population stratification is a major source of inflated false positive rates in genome wide association studies. However, relatively few studies have examined its impact on gene-gene interaction detection, despite the importance of epistasis for understanding the genetic architecture of complex traits. In this study, we identify scenarios under which population stratification can inflate the interaction test statistics. Through analytical derivations and simulation studies, we show that this inflation is not adequately controlled by including principal components as covariates in the regression model. We then propose an alternative approach that effectively controls the inflation of false-positive rates for interaction test statistics due to population stratification by using single nucleotide polymorphism-by-population structure interaction as an additional covariate term in the regression model.

Predicting implantation outcomes in vitro fertilization (IVF) is challenging because clinically relevant information becomes available at different points during treatment. In this study, we developed a stage-aware, leakage-conscious, and explainable machine-learning framework to predict implantation using a public IVF dataset. {beta}-HCG was excluded to avoid potential downstream target leakage. The data were divided into training, validation, and independent test sets. Several machine-learning models, including Logistic Regression, Random Forest (RF), Extra Trees, XGBoost, and Hist Gradient Boosting, were optimized using stratified cross-validation. A fixed RF model was then evaluated at three stages of the IVF process, using baseline couple characteristics, gamete and fertilization variables, and embryological variables, respectively. Model performance improved as additional treatment-related information became available, with the largest improvement observed after the inclusion of embryological variables. The final RF model showed good predictive performance on the independent test set. The number of embryos was consistently identified as the most important predictor across RF, permutation, and SHAP analyses. Overall, the findings suggest that IVF implantation prediction becomes more accurate as later-stage treatment information is incorporated, particularly information related to embryological characteristics. The proposed approach also provides interpretable predictions while reducing the risk of information leakage. Further evaluation using independent clinical datasets is needed before the framework can be considered for clinical use.

In competing risks settings, covariate effects and group comparisons are usually assessed one event at a time - through log-rank or Cox tests on the cause-specific hazards, or Gray's test or Fine-Gray regression on a cumulative incidence function (CIF). This can obscure a clinically important quantity: the ratio between the event of interest and the competing event, since groups may differ little on the individual events yet differ sharply in their ratio. The generalized competing event (GCE) framework makes this ratio the object of inference; on the cause-specific scale the hazard ratio omega+(t) = lambda_1(t)/lambda_2(t) is estimated efficiently from a single stacked (Lunn-McNeil) model. We extend the framework to two scales that describe realized incidence. The subdistribution hazard ratio omega-tilde+(t) = lambda-tilde_1(t)/lambda-tilde_2(t) is estimated by a stacked, risk-set-weighted extension of the Lunn-McNeil construction; the cumulative-incidence ratio rho(t) = F_1(t)/F_2(t) - the odds that a subject's realized event by time t is the event of interest - by jackknife pseudo-observation regression of the Aalen-Johansen estimator. We relate the three contrasts: rho equals omega+ exactly under proportional cause-specific hazards, and equals omega-tilde+ only in the small-time limit under proportional subdistribution hazards, drifting toward 1 thereafter. The orthogonality that makes omega+ efficient is lost on both cumulative-incidence scales - omega tilde+ through overlapping weighted risk sets and shared censoring weights, rho through the shared all-cause survivor - so each carries a covariance term that must be handled and that bounds efficiency relative to the hazard-scale test. We derive the corresponding variances, study operating characteristics by simulation, illustrate on hypothetical prostate and head-and-neck cohorts, and provide an implementation in the gcemod R package.

Dysfunctions in transcription cause severe human diseases. In eukaryotes, transcription initiation by RNA polymerase (Pol) II requires the formation of the pre-initiation complex (PIC) composed of general transcription factors (GTFs), Pol II, and an essential multiprotein coactivator, Mediator. In vitro experiments with GTFs and Pol II have depicted a linear sequence of PIC assembly. However, PIC assembly paths in vivo at the genomic scale and Mediator's role in these paths remain poorly understood. One strategy characterizes ChIP-seq occupancy alterations in Mediator mutants, but analysis of these data still lacks a formalized interpretation framework. In this study, we address this need by developing a quantitative framework that scores candidate assembly paths according to their ability to explain these alterations. This framework combines a model of PIC assembly based on a system of differential equations with parameter optimization based on a penalized criterion. We applied this new framework to datasets from two yeast Mediator mutants, in the Med17 and Med10 subunits. Our results show that the best-supported PIC assembly paths differ from those inferred in vitro. In particular, we show that it is important to distinguish the core TFIIH module from the TFIIK kinase module, and that similar paths starting with Mediator or TBP receive different levels of support. Furthermore, although multiple paths fit the data equally well, they differ markedly for specific gene groups. Taken together, the results show that the path can depend on the gene, and are also compatible with more than one path being operational in vivo for a single gene. This work improves our understanding of PIC assembly paths and transcription initiation, and opens perspectives for future modeling and experimental studies.

Understanding the relationship between gene expression dynamics and cellular identity remains a central challenge in single cell biology. Here, we introduce a novel computational and mathematical framework that integrates information geometry, fuzzy topology, and UMAP analysis to model gene expression landscapes derived from single cell RNA sequencing data. We formalize gene expression data as a fuzzy topological space, where interactions between expression points are governed by probabilistic distributions inspired by manifold learning approaches such as UMAP. Within this framework, we define an information geometric structure through a Fisher metric induced by these distributions, enabling the computation of geodesic trajectories that capture cellular differentiation processes. A key contribution of this work is the derivation of analytical conditions, expressed as expression radius formulas, that characterize local neighborhoods in gene expression space. These conditions allow for the identification of genes associated with stem cell states and predictions in transitional cell types in future work. Application of the proposed framework to single cell datasets reveals biologically meaningful gene sets enriched in key regulatory pathways and transcription factors, demonstrating the capacity of our approach to uncover latent structure in complex gene expression data. Our results suggest that integrating differential geometry with statistical learning theory offers a powerful paradigm for modeling genotype and phenotype relationships and cellular state transitions, with potential implications for precision medicine and systems biology.

Diffusion within the cytoplasm is fundamental to numerous biological processes. Fluorescence recovery after photobleaching (FRAP) is one of the most common method for quantifying molecular diffusivity in living cells using standard laser scanning confocal microscopy (LSCM). However, accurately measuring fast cytoplasmic diffusion (typically >10 m^2/s) is challenging due to rapid recovery kinetics, weak signal-to-noise ratios, post-bleach signal artifacts, and spatial restrictions affecting normalization. While individual challenges have been addressed in specific contexts, a simple and robust framework to quantify cytoplasmic diffusivity remains elusive. Here, we present a FRAP methodology specifically designed to overcome these obstacles. By utilizing the Gaussian function -- the impulse response (ImpRes) of the diffusion equation in an infinite medium -- our approach leverages the full spatiotemporal dataset through a single-equation three-parameter fitting procedure, thus releasing restrictions to small regions of interest and arbitrary initial time-points. The methodology was validated on three datasets of increasing complexity: in silico simulated recovery profiles, in vitro data from FITC-dextran in glycerol solution, and live-cell imaging of free cytoplasmic GFP. Systematic comparison with existing models demonstrates that the ImpRes approach significantly reduces sensitivity to noise and imperfect fluorescence normalization, while remaining robust against short-term biases, such as transient probe photo-activation. Given its robustness under realistic experimental conditions and its ease of implementation, the proposed FRAP methodology provides a reliable tool for quantitative cytoplasmic analysis.

Passive acoustic sensors and machine learning classifiers offer a scalable approach to monitoring wildlife populations. These classifiers typically detect whether a species is present in a short time-window of audio; ecologists, however, want ecologically meaningful measures such as indices of site-level abundance. Vocal density, the proportion of time-windows containing a vocalization, is a useful abundance proxy. Recovering vocal density from classifiers requires calibrating their outputs against expert-labeled data. Because classifier performance varies between sites, a single global calibration yields overconfident site-level estimates, yet calibrating every site independently requires prohibitive labeling effort. Estimating site-level vocal density at the scale of modern sensor networks demands a label-efficient alternative. We present a Bayesian hierarchical extension of Platt scaling, a commonly used calibration method, that calibrates classifiers at the site level while borrowing strength across sites. Where site-level covariates are available, a Gaussian-process prior lets the model learn how calibration varies across covariate space. Averaging calibrated per-clip probabilities across a site's recordings yields a posterior distribution over vocal density, and we show how to carry that uncertainty into downstream regressions of vocal density on environmental covariates. Across simulations and two fully annotated field datasets from Hawai'i and the Pacific Northwest, our model attained desired credible-interval coverage of vocal density estimates with narrower intervals than independent site-level calibration at equal labeling effort, and lower mean squared error. Global calibration, by contrast, failed to attain coverage except when simulated sites were genuinely homogeneous. Unlike site-level calibration, our model enables calibration at sites with zero labeled data, and incorporating covariates further improved precision when site heterogeneity was covariate-driven. Applied to a 283-site dataset in Pennsylvania with only two labeled clips per site, we recovered known habitat associations for the declining Wood Thrush (Hylocichla mustelina). Our method provides a principled, label-efficient route from bioacoustic classifier scores to site-level abundance indices with well-quantified uncertainty, making rigorous ecological inference feasible at the scale at which acoustic sensor networks are now deployed.

Accurate normalization is essential for differential expression analysis of RNA-sequencing data. Popular normalization methods such as the median-of-ratios and trimmed mean of M-values do not leverage information from the experimental design. This may be inefficient in experiments with large-scale systematic expression changes or complex designs. Here, we introduce design-informed size factor estimation (disize), a normalization method that uses information from the experimental design to improve accuracy. disize uses a modified generalized linear mixed model to robustly distinguish between biological signal and sample-specific size factors. We also propose a mechanistically justified data-generating process for RNA-sequencing counts that is derived from previous models of transcription and sequencing. Through simulations based on this data-generating process and validating on true RNA-seq data, we show that disize recovers size factors more accurately than existing methods, particularly in challenging scenarios with low gene expression and a high proportion of differentially expressed genes; this in turn improves downstream analysis. disize provides a robust and accurate approach to normalization, highlighting the significant benefits of integrating experimental design information directly into normalization for transcriptomic datasets. Author summaryIn transcriptomic analysis, normalization adjusts for technical biases arising from library preparation and sequencing. Methods implemented in widely used packages like DESeq2 and edgeR ignore information in the experimental design during normalization. Incorporating information from the experimental design into a normalization method has the potential to yield more accurate results. To do this, we developed a new method, design-informed size factor estimation (disize), that uses a statistical model to jointly account for the biological signal defined by the design and the sample-specific batch effect. By separating the biological variation into its components, disize can more robustly estimate the batch effect. To validate our approach, we constructed a flexible simulation framework relying on a mechanistically justified data-generating process for RNA-seq data. Our benchmarks on both simulated and true RNA-seq data show that disize recovers the true size factors more accurately than existing methods, particularly in challenging scenarios with low counts or a high proportion of differentially expressed genes. This improved normalization yields more reliable downstream results in differential expression analysis.